!
!  ifits    I/O FITS interface 
!
!  Copyright © 2010 - 2011 F.Hroch (hroch@physics.muni.cz)
!
!  This file is part of Munipack.
!
!  Munipack is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!  
!  Munipack is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  
!  You should have received a copy of the GNU General Public License
!  along with Munipack.  If not, see <http://www.gnu.org/licenses/>.


module ifits

  use fitsio

  integer, parameter, private :: rp = selected_real_kind(15)

contains

subroutine finit(filename,append,raper,ring)

  implicit none

  character(len=*),intent(in) :: filename
  logical, intent(in) :: append
  real, dimension(:),intent(in) :: raper,ring
  integer :: status, blocksize, hdutype, extver = 0
  character(len=80) :: buf, label, label_err
  character(len=80), dimension(11+2*size(raper)) :: ttype, tform, tunit
  integer :: i,j

  ttype(1) = 'x';            tform(1) = '1D';  tunit(1) = ''
  ttype(2) = 'y';            tform(2) = '1D';  tunit(2) = ''
  ttype(3) = 'sky';          tform(3) = '1D';  tunit(3) = ''
  ttype(4) = 'sky_err';      tform(4) = '1D';  tunit(4) = ''
  ttype(5) = 'peak_ratio';   tform(5) = '1D';  tunit(5) = ''
  ttype(6) = 'apcflux';      tform(6) = '1D';  tunit(6) = ''
  ttype(7) = 'apcflux_err';  tform(7) = '1D';  tunit(7) = ''
  ttype(8) = 'psfflux';      tform(8) = '1D';  tunit(8) = ''
  ttype(9) = 'psfflux_err';  tform(9) = '1D';  tunit(9) = ''
  ttype(10) = 'psfpeak';     tform(10) = '1D'; tunit(10) = ''
  ttype(11) = 'psfpeak_err'; tform(11) = '1D'; tunit(11) = ''

  do i = 1, size(raper)

     write(buf,'(i0)') i
     label = 'apflux'//trim(buf)
     label_err= trim(label)//'_err'
     j = 10 + 2*i
!     write(*,*) trim(label),' ',trim(label_err)
     ttype(j)   = trim(label);     tform(j) = '1D';   tunit(j) = ''
     ttype(j+1) = trim(label_err); tform(j+1) = '1D'; tunit(j+1) = ''

  end do


  status = 0
  if( append ) then
     call ftopen(15,filename,1,blocksize,status)
  else
     blocksize = 0
     call ftinit(15,filename,blocksize,status)
  end if

  if( status /= 0 ) then
     write(*,*) 'Failed to init data store file: `',trim(filename),"'"
     stop 22
  end if

  ! look for the extension
  call ftmnhd(15,BINARY_TBL,MEXTNAME,extver,status)

  if( status == 0 ) then
     call ftdhdu(15,hdutype,status)
     if( status /= 0 ) write(*,*) 'Failed to delete MUNIPACK extension.'
  end if

  status = 0

  call ftibin(15,0,size(ttype),ttype,tform,tunit,MEXTNAME,0,status)

  do i = 1, size(raper)
     write(label,'(a,i0)') 'APER',i
     call ftpkye(15,trim(label),raper(i),-5,'[pix] aperture radius',status)
  end do

  call ftpkye(15,'ANNULUS1',ring(1),-5,'[pix] inner sky annulus radius',status)
  call ftpkye(15,'ANNULUS2',ring(2),-5,'[pix] outer sky annulus radius',status)
  call ftclos(15,status)

  if( status /= 0 )then
     write(*,*) 'FITSIO: following error(s) has been occured:'
     do
        Call ftgmsg(buf)
        if( buf == ' ' ) exit
        write(*,*) trim(buf)
     enddo
  endif


end subroutine finit

subroutine ffind(filename,output,maxsky,saturation,fwhm,lothresh, &
     threshold,shrplo,shrphi,rndlo,rndhi,readns, phpadu)

!  use cfitsio
!  use dao
  use mdaofind

  implicit none

  integer, parameter :: DIM = 2

  character(len=*),intent(in) :: filename, output
  real,intent(in) :: saturation,fwhm,lothresh,threshold,&
       shrplo,shrphi,rndlo,rndhi,readns, phpadu
  integer,intent(in) :: maxsky
  integer :: status, blocksize, naxis, pcount,gcount, bitpix
  integer :: group = 1, extver = 0, frow = 1, felem = 1
  real :: nullval = 0.0
  integer, dimension(DIM) :: naxes
  logical :: simple,extend,anyf
  real, dimension(:,:), allocatable :: data
  real, dimension(:), allocatable :: xstar,ystar,hstar!,xs,ys,hs
!  integer, dimension(:), allocatable :: index
  character(len=80) :: buf
  real :: lobad, aduthresh
  logical :: alloc = .false.
!  integer :: i,idx
!  character(len=80), dimension(12) :: ttype, tform, tunit

!  ttype(1) = 'x';            tform(1) = '1D';  tunit(1) = ''
!  ttype(2) = 'y';            tform(2) = '1D';  tunit(2) = ''
!  ttype(3) = 'apflux';       tform(3) = '1D';  tunit(3) = ''
!  ttype(4) = 'apflux_err';   tform(4) = '1D';  tunit(4) = ''
!  ttype(5) = 'apcflux';      tform(5) = '1D';  tunit(5) = ''
!  ttype(6) = 'apcflux_err';  tform(6) = '1D';  tunit(6) = ''
!  ttype(7) = 'psfflux';      tform(7) = '1D';  tunit(7) = ''
!  ttype(8) = 'psfflux_err';  tform(8) = '1D';  tunit(8) = ''
!  ttype(9) = 'psfpeak';      tform(9) = '1D';  tunit(9) = ''
!  ttype(10) = 'psfpeak_err'; tform(10) = '1D'; tunit(10) = ''
!  ttype(11) = 'sky';         tform(11) = '1D'; tunit(11) = ''
!  ttype(12) = 'sky_err';     tform(12) = '1D'; tunit(12) = ''
!
!  write(*,*) trim(filename)
  
  status = 0
  call ftopen(15,filename,0,blocksize,status)
  call ftghpr(15,DIM, &
       simple,bitpix,naxis,naxes,pcount,gcount,extend,status)
  allocate(data(naxes(1),naxes(2)))
  call ftg2de(15,group,nullval,size(data,1),naxes(1),naxes(2), &
       data,anyf,status)
  call ftclos(15,status)
  if( status /= 0 ) stop 'Failed to read image data.'

     !  write(*,*) data(23,78)

  
  call ftopen(15,output,1,blocksize,status)

  call ftmnhd(15,BINARY_TBL,MEXTNAME,extver,status)
  if( status == BAD_HDU_NUM ) &
       write(*,*) 'MUNIPACK extension not found in find. Uninitialized?'
!       stop 'MUNIPACK extension not found in find. Uninitialized?'

!  call finder(data,threshold,fwhm,maxsky,saturation,xstar,ystar)
!  write(*,*) size(data,1),size(data,2)

  if( status == 0 ) then
     
     alloc = .true.

     call daofind(data,maxsky,saturation,fwhm,lothresh,threshold, &
          shrplo,shrphi,rndlo,rndhi,readns,phpadu,xstar,ystar,hstar, &
          lobad, aduthresh)

     ! sort arrays by height above lower threshold
     call sorter(xstar,ystar,hstar)

  end if

!  allocate(index(size(hstar)))
!  forall( i = 1:size(hstar) )
!     index(i) = i
!  end forall
!  call quicksorti(hstar,size(hstar),index)
!  allocate(xs(size(xstar)),ys(size(ystar)),hs(size(hstar)))
!  xs = xstar
!  ys = ystar
!  hs = hstar
!  do i = 1, size(index)
!     idx = index(size(index) - i + 1)
!     xstar(i) = xs(idx)
!     ystar(i) = ys(idx)
!  end do
!  do i = 1,size(hstar)/2
!     hstar(i) = hs(size(index) - i + 1)
!  end do
!  deallocate(index,xs,ys,hs)
!
!  write(*,*) hstar
!  write(*,*) index


! results fills FITS extension
  call ftpcom(15,'Parameters used during star detection:',status)

  write(buf,*) 'Pixels used for sky determination =',maxsky
  call ftpcom(15,buf,status)

!  write(buf,*) 'Saturation =',saturation
!  call ftpcom(15,buf,status)

  write(buf,*) 'FWHM (pix) =',fwhm
  call ftpcom(15,buf,status)

  write(buf,*) 'Threshold (sigma)=',threshold
  call ftpcom(15,buf,status)

  write(buf,*) 'Lower threshold (sigma)=',lothresh
  call ftpcom(15,buf,status)

!  write(buf,*) ' =',lothresh
!  call ftpcom(15,buf,status)

! uchovat tyhle parametry:
!  call wrhead (3, 1, ncol, nrow, 6, lobad, hibad, hmin, 0., phpadu, readns, 0.)


!  write(buf,*) 'Parameters: threshold=',threshold,' fwhm =',fwhm
!  call FTPCOM(15,buf,status)

  call ftukye(15,'GAIN',phpadu,-5,'[e-/ADU] gain',status)
  call ftukye(15,'LOWBAD',lobad,-5,'[ADU] low good datum',status)
  call ftukye(15,'HIGHBAD',saturation,-5,'[ADU] high good datum',status)
  call ftukye(15,'RNOISE',readns,-5,'[ADU] read noise',status)
  call ftukye(15,'THRESH',aduthresh,-5,'[ADU] threshold',status)

!  write(*,*) size(xstar),xstar

  call ftpcle(15,1,frow,felem,size(xstar),xstar,status)
  call ftpcle(15,2,frow,felem,size(ystar),ystar,status)
  call ftpcle(15,5,frow,felem,size(hstar),hstar,status)

  call ftclos(15,status)

  if( alloc ) deallocate(data,xstar,ystar,hstar)

  if( status /= 0 )then
     write(*,*) 'FITSIO: following error(s) has been occured:'
     do
        Call ftgmsg(buf)
        if( buf == ' ' ) exit
        write(*,*) trim(buf)
     enddo
  endif


contains

  subroutine sorter(xstar,ystar,hstar)

    use quicksort

    real, dimension(:),intent(inout) :: xstar,ystar,hstar

    real(rp), dimension(:), allocatable :: xs,ys,hs
    integer, dimension(:), allocatable :: id
    integer :: i,n,idx

    allocate(xs(size(xstar)),ys(size(ystar)),hs(size(hstar)),id(size(hstar)))
    forall( i = 1:size(hstar) )
       id(i) = i
       hs(i) = hstar(i)
    end forall
    call qsort(hs,id)
    ! sorted from low to high

    xs = xstar
    ys = ystar
    n = size(id) + 1
    do i = 1, size(id)
       idx = id(n - i)
       xstar(i) = xs(idx)
       ystar(i) = ys(idx)
    end do

    ! reverse sort
    forall( i = 1:size(hstar) )
       hstar(i) = hs(n - i)
    end forall
    deallocate(id,xs,ys,hs)
    
  end subroutine sorter


end subroutine ffind


subroutine faphot(filename,output)!,raper,ring)

!  use cfitsio
!  use dao
  use mdaofotometr

  implicit none

  integer, parameter :: DIM = 2

  character(len=*),intent(in) :: filename,output
  integer :: status, blocksize, naxis, pcount,gcount, bitpix, nrows
  integer :: group = 1, extver = 0, frow = 1, felem = 1
  real :: nullval = 0.0
  integer, dimension(DIM) :: naxes
  logical :: simple,extend,anyf
  real, dimension(:,:), allocatable :: data
  real, dimension(:), allocatable :: xstar,ystar,hstar,apcflux,apcflux_err, &
       sky,sky_err, apf
  real, dimension(:,:), allocatable :: apflux,apflux_err
  real, dimension(2) :: ring
  real, dimension(12) :: raper
  real :: lobad, hibad, phpadu
  logical :: alloc = .false.
  integer :: i
  character(len=80) :: buf

  status = 0
  call ftopen(15,filename,0,blocksize,status)
  call ftghpr(15,DIM, &
       simple,bitpix,naxis,naxes,pcount,gcount,extend,status)
  allocate(data(naxes(1),naxes(2)))
  call ftg2de(15,group,nullval,size(data,1),naxes(1),naxes(2), &
       data,anyf,status)
!  if( status /= 0 ) stop 'Failed to read image data.'
  call ftclos(15,status)

  call ftopen(15,output,1,blocksize,status)
  call ftmnhd(15,BINARY_TBL,MEXTNAME,extver,status)
  if( status == BAD_HDU_NUM ) &
       stop 'MUNIPACK extension not found in aphot. Uninitialized?'

  call ftgnrw(15,nrows,status)

  allocate(xstar(nrows),ystar(nrows),hstar(nrows))
  allocate(apflux(nrows,size(raper)),apflux_err(nrows,size(raper)))
  allocate(apcflux(nrows),apcflux_err(nrows),sky(nrows),sky_err(nrows))
  
  call ftgcve(15,1,frow,felem,size(xstar),nullval,xstar,anyf,status)
  call ftgcve(15,2,frow,felem,size(ystar),nullval,ystar,anyf,status)
  call ftgcve(15,5,frow,felem,size(hstar),nullval,hstar,anyf,status)

  call ftgkye(15,'LOWBAD',lobad,buf,status)
  if( status /= 0 ) lobad = 0.0
  status = 0
  call ftgkye(15,'HIGHBAD',hibad,buf,status)
  if( status /= 0 ) hibad = huge(0.0)
  status = 0
  call ftgkye(15,'GAIN',phpadu,buf,status)
  if( status /= 0 ) phpadu = 1.0
  status = 0

  ring = (/20.0, 30.0/)
  call ftgkye(15,'ANNULUS1',ring(1),buf,status)
  status = 0
  call ftgkye(15,'ANNULUS2',ring(2),buf,status)
  status = 0

  do i = 1, size(raper)
     write(buf,'(a,i0)') 'APER',i
     call ftgkye(15,buf,raper(i),buf,status)
     status = 0
  end do

  apflux = -1.0
  apflux_err = -1.0

  if( status == 0 ) then

     alloc = .true.

     call daophotsb(data,xstar,ystar,hstar,raper,ring,lobad,hibad,phpadu, &
          apflux,apflux_err,apcflux,apcflux_err,sky,sky_err)
  end if

!  call ftpcom(15,'Parameters used during aperture photometry:',status)

!  write(buf,*) 'Ring radii: inner=',ring(1),' outer=',ring(2)
!  call ftpcom(15,buf,status)

!  do i = 1, size(raper)

!     write(buf,*) 'Aperture radius ',i,'= ',raper(i)
!     call ftpcom(15,buf,status)

!  end do

  do i = 1,size(apflux,2)
     allocate(apf(size(apflux(:,i))))
     apf = apflux(:,i)
     call ftpcle(15,10+2*i,frow,felem,size(apf),apf,status)
     apf = apflux_err(:,i)
     call ftpcle(15,10+2*i+1,frow,felem,size(apf),apf,status)
     deallocate(apf)
  end do

  call ftpcle(15,3,frow,felem,size(sky),sky,status)
  call ftpcle(15,4,frow,felem,size(sky_err),sky_err,status)

  call ftclos(15,status)

  if( alloc ) &
       deallocate(data,xstar,ystar,apflux,apflux_err,apcflux,apcflux_err, &
       sky,sky_err)

  if( status /= 0 )then
     write(*,*) 'FITSIO: following error(s) has been occured:'
     do
        Call ftgmsg(buf)
        if( buf == ' ' ) exit
        write(*,*) trim(buf)
     enddo
  endif


end subroutine faphot

end module ifits
